Para esta práctica he implementado un workflow en el que encontramos diferentes modulos:
Volcano plot
En este trabajo explicaré los diferentes resultados así como las decisiones que se han obtenido a lo largo del proyecto. La idea es generar un workflow que nos devuelva información útil para tratar datos de BRCA obtenidos de TCGA. Concretamente, busco que el flujo de trabajo determine si merece la pera utilizar Lasso/ridge o algún otro método de machine learning. Además vamos a buscar conexiones funcionales que podrían dar sentido biológico a los genes que hay en cada modelo predictivo a través de un enriquecimiento funcional de dichos genes. Y finalmente comprobar si podemos agrupar los diferentes pacientes en diferentes clusters y analizar la diferencia de supervivencia de los diferentes subgrupos.
El workflow ha sido implementado en bash. Podemos seleccionar diferentes parámetros a la hora de ejecutarlo. Contiene una función usage en caso de que se ecriban mal los parámetros para orientarnos sobre como deberíamos escribirlos.
usage: ./work.sh [-svb] [-l number] [-p number] [-t number]
-s skip downloading the data if downloaded
-v create volcano plot
-b Balancea los datos para las predicciones
-l number Set LogFC, for downloading the data, for the volcano or both.
-f string Set the feature selection method. (NONE, COX or RFI)
-p number Set P-value, for downloading the data, for the volcano or both.
-t number Set the step for the data replication.
Los resultados que me han parecido interesantes guardar se encuentrar organizados dentro de la carpeta results mientras que los datos intermedios y datos brutos se encuentran en la carpeta de datos.
Para seleccionar las variables de los datos clínicos, he comprobado que no tengan muchos valores perdidos y he buscado aquellas variables que resulten mas informativas. He decidido truncar los pacientes con valores perdidos; como posteriormente los datos se van a replicar podemos permitirnos truncar los pacientes que tengar valores perdidos.
Concretamente se han seleccionado las siguientes variables:
Estos datos clínicos, como podemos apreciar, han sido obtenidos de una tabla donde encotramos muestras de cancer de mama clasificadas por subtipos de PAM50, y por ER status, PR status y HER2 Final status.
Los datos de expresión han sido obtenidos de TCGA como vimos en clase:
brcaData <- getFirehoseData(dataset="BRCA", runDate="20160128",gistic2Date="20160128",forceDownload=F, clinical =TRUE, RNASeq2GeneNorm =TRUE)
brca_rnaseq <- getData(brcaData,type = "RNASeq2GeneNorm")
brca_rnaseq.tumour <- brca_rnaseq[, which(as.numeric(substr(colnames(brca_rnaseq), 14,15)) < 10)]
colnames(brca_rnaseq.tumour) <- substr(colnames(brca_rnaseq.tumour), 1,12)
brca_rnaseq.tumour <- brca_rnaseq.tumour[, !duplicated(colnames(brca_rnaseq.tumour))]
Una vez descargados hemos seleccionados aquellos genes que aparecen tanto en Luminal como en TNBC. El subtipo de Luminal tiene baja proliferación y tiene menos recaida, mientras que el subtipo de TNBC es mucho mas agresivo. Los resultados de expresión diferencial para cada subtipo se van a unir generando una matriz de diseño.
A continuación haremos el análisis de expresión diferencial. Convertimos nuestra matriz de counts para poder aplicar el método de Voom que estima la variación media de los log-counts (de counts) generando pesos para cada observación. Así ya tendremos los datos preparados para poder aplicar a estos datos de expresión los métodos ideados para el análisis de microarrays.
tnbc_samples <- sample_data %>% dplyr::filter(`ER Status` == "Negative" & `PR Status` == "Negative" & `HER2 Final Status` == "Negative" & `PAM50 mRNA` != "Luminal A")
tnbc_barcodes <- tnbc_samples$`Complete TCGA ID`
brca_rnaseq.tnbc <- brca_rnaseq.tumour[, which(colnames(brca_rnaseq.tumour) %in% tnbc_barcodes)]
luminal_samples <- sample_data %>% dplyr::filter(`PAM50 mRNA` == "Luminal A")
luminal_barcodes <- luminal_samples$`Complete TCGA ID`
brca_rnaseq.luminal <- brca_rnaseq.tumour[, which(colnames(brca_rnaseq.tumour) %in% luminal_barcodes)]
rnaseq.for.de <- cbind(brca_rnaseq.luminal, brca_rnaseq.tnbc)
counts = rnaseq.for.de[apply(rnaseq.for.de,1,function(x) sum(x==0))<ncol(rnaseq.for.de)*0.8,]
df.l <- data_frame("sample" = colnames(brca_rnaseq.luminal), "status" = rep(0, length(colnames(brca_rnaseq.luminal))) )
df.t <- data_frame("sample" = colnames(brca_rnaseq.tnbc), "status" = rep(1, length(colnames(brca_rnaseq.tnbc))) )
df <- rbind(df.l,df.t)
design <- model.matrix(~ status, data = df)
dge <- DGEList(counts=counts)
A <- rowSums(dge$counts)
isexpr <- A > 100 # Keeping genes with total counts more than 100.
dge <- calcNormFactors(dge)
v <- voom(dge[isexpr,], design, plot=FALSE)
fit <- lmFit(v, design)
fit <- eBayes(fit)
diff.exp.df <- topTable(fit, coef = "status", n = Inf, sort = "p", p = 0.01) # Positive log-fold-changes mean higher expression in d1
diff.exp.df$gene.name <- rownames(diff.exp.df)
Ahora ya podemos sacar los valores de la expresión diferencial de los genes. Una vez tenemos todos los datos, mostraremos los siguientes sets de datos.
De esta forma podemos tener un seguimiento rápido sobre el filtrado general que hemos hecho a los datos.
menor <- diff.exp.df[diff.exp.df$logFC<(-lfc),]
mayor <- diff.exp.df[diff.exp.df$logFC>lfc,]
sobreexp <- rbind(menor,mayor)
sobreexp <- sobreexp[ sobreexp$adj.P.Val<pval,]
genes.model <- counts[sobreexp$gene.name,]
genes.model <- t(genes.model)
genes.model<- data.frame(genes.model)
genes.model$`Complete TCGA ID` <- rownames(genes.model)
compact<- sample_data[c(sample_data$`Complete TCGA ID`) %in% genes.model$`Complete TCGA ID`,]
join <- inner_join(compact, genes.model)
clinicos.sin.quitar.vacios<- sample_data_1
clinicos.completos.sin.exp <- sample_data
clinicos.reducidos <- compact
clinicos.expresion <- join
cat("Los datos originales tienen unas dimensiones de: ", dim(clinicos.sin.quitar.vacios)[1],"x",dim(clinicos.sin.quitar.vacios)[2],
"\nLos datos filtrados tienen unas dimensiones de: ", dim(clinicos.completos.sin.exp)[1],"x",dim(clinicos.completos.sin.exp)[2],
"\nLos datos reducidos tienen unas dimensiones de: ", dim(clinicos.reducidos)[1],"x",dim(clinicos.reducidos)[2],
"\nLos datos finales tienen unas dimensiones de: ", dim(clinicos.expresion)[1],"x",dim(clinicos.expresion)[2], "\n"
)
## Los datos originales tienen unas dimensiones de: 825 x 30
## Los datos filtrados tienen unas dimensiones de: 515 x 15
## Los datos reducidos tienen unas dimensiones de: 311 x 15
## Los datos finales tienen unas dimensiones de: 311 x 1043
Los datos que utiliza son los de expresión de los genes diff.exp.df. Para hacer el volcano plot sencillamente cogeremos los genes que aparecen expresados diferencialmente y los plotearemos como se ve en el codigo.
cat("El lfc para los datos descargados es ", lfc," y el p-valor es de ", pval, "\n")
## El lfc para los datos descargados es 4 y el p-valor es de 0.01
tab = data.frame(logFC = diff.exp.df$logFC, negLogPval= -log10(diff.exp.df$adj.P.Val))
tab2 = data.frame(logFC = diff.exp.df$logFC, negLogPval= -log10(diff.exp.df$adj.P.Val), Gene=diff.exp.df$gene.name)
par(mar = c(5, 4, 4, 5))
plot(tab, pch = 16, cex = 0.6, xlab = expression(log[2]~fold~change), ylab = expression(-log[10]~pvalue))
#signGenes = (abs(tab$logFC) > lfc & tab$negLogPval > -log10(pval))
points(tab[(abs(tab$logFC) > lfc), ], pch = 16, cex = 0.8, col = "orange")
points(tab[(tab$negLogPval > -log10(pval)), ], pch = 16, cex = 0.8, col = "green")
points(tab[(abs(tab$logFC) > lfc & tab$negLogPval > -log10(pval)), ], pch = 16, cex = 0.8, col = "red")
abline(h = -log10(pval), col = "green3", lty = 2)
abline(v = c(-lfc, lfc), col = "blue", lty = 2)
mtext(paste("pval =", pval), side = 4, at = -log10(pval), cex = 0.8, line = 0.5, las = 1)
mtext(c(paste("-", lfc, "fold"), paste("+", lfc, "fold")), side = 3, at = c(-lfc, lfc), cex = 0.8, line = 0.5)
with(subset(tab2, negLogPval > -log10(pval) & abs(logFC)>lfc), textxy(logFC, negLogPval, labs=Gene, cex=.4))
Como vemos este gráfico resulta muy ilustrativo para poder apreciar cuantitativamente la expresión diferencial de los genes.
El siguiente modulo que nos encontramos en el workflow es el de selección de variables y enriquecimiento de genes. Lo se ha hecho en este módulo es consultar en la ontología de funciones por el conjunto de genes que tenemos. El codigo que tenemos a continuación se usará dos veces, la primera vez nos devolverá el enriquecimiento de los genes que tenemos expresados diferencialmente.
enrich_cp = function(res, comparison, type="over") {
res = res %>% data.frame() %>% left_join(entrezsymbol, by = "hgnc_symbol") %>% filter(!is.na(entrezgene))
# universe = brcaData@GISTIC@AllByGene$Gene.Symbol
if(type=="all"){
res <- res %>% filter(abs(logFC) > lfc & adj.P.Val < pval) # lfc and pval threshold defined above in the volcano plot
genes = res$entrezgene
mf = enrichGO(genes, OrgDb = orgdb, ont = "MF", pAdjustMethod = "BH",
qvalueCutoff = 1, pvalueCutoff = 1)
cc = enrichGO(genes, OrgDb = orgdb, ont = "CC", pAdjustMethod = "BH",
qvalueCutoff = 1, pvalueCutoff = 1)
bp = enrichGO(genes, OrgDb = orgdb, ont = "BP", pAdjustMethod = "BH",
qvalueCutoff = 1, pvalueCutoff = 1)
kg = enrichKEGG(gene = genes, organism = keggname, pvalueCutoff = 1,
qvalueCutoff = 1, pAdjustMethod = "BH")
all = list(mf = mf, cc = cc, bp = bp, kg = kg)
all[["summary"]] = summarize_cp(all, comparison)
return(all)
}
if(type=="over"){
res.over <- res %>% filter(logFC > lfc & adj.P.Val < pval)
genes = res.over$entrezgene
mf = enrichGO(genes, OrgDb = orgdb, ont = "MF", pAdjustMethod = "BH",
qvalueCutoff = 1, pvalueCutoff = 1)
cc = enrichGO(genes, OrgDb = orgdb, ont = "CC", pAdjustMethod = "BH",
qvalueCutoff = 1, pvalueCutoff = 1)
bp = enrichGO(genes, OrgDb = orgdb, ont = "BP", pAdjustMethod = "BH",
qvalueCutoff = 1, pvalueCutoff = 1)
kg = enrichKEGG(gene = genes, organism = keggname, pvalueCutoff = 1,
qvalueCutoff = 1, pAdjustMethod = "BH")
all = list(mf = mf, cc = cc, bp = bp, kg = kg)
all[["summary"]] = summarize_cp(all, comparison)
return(all)
}
if(type=="under"){
res.under <- res %>% filter(logFC < -lfc & adj.P.Val < pval)
genes = res.under$entrezgene
mf = enrichGO(genes, OrgDb = orgdb, ont = "MF", pAdjustMethod = "BH",
qvalueCutoff = 1, pvalueCutoff = 1)
cc = enrichGO(genes, OrgDb = orgdb, ont = "CC", pAdjustMethod = "BH",
qvalueCutoff = 1, pvalueCutoff = 1)
bp = enrichGO(genes, OrgDb = orgdb, ont = "BP", pAdjustMethod = "BH",
qvalueCutoff = 1, pvalueCutoff = 1)
kg = enrichKEGG(gene = genes, organism = keggname, pvalueCutoff = 1,
qvalueCutoff = 1, pAdjustMethod = "BH")
all = list(mf = mf, cc = cc, bp = bp, kg = kg)
all[["summary"]] = summarize_cp(all, comparison)
return(all)
}
}
convert_enriched_ids = function(res, entrezsymbol) {
res = res %>% mutate(geneID = strsplit(as.character(geneID), "/")) %>% tidyr::unnest(geneID) %>%
left_join(entrezsymbol, by = c(geneID = "entrezgene")) %>% group_by(ID,
Description, GeneRatio, BgRatio, pvalue, p.adjust, qvalue, Count, ont,
comparison) %>% summarise(geneID = paste(geneID, collapse = "/"), symbol = paste(hgnc_symbol,
collapse = "/"))
return(res)
}
orgdb = "org.Hs.eg.db"
biomart_dataset = "hsapiens_gene_ensembl"
keggname = "hsa"
mart = biomaRt::useMart(biomart = "ensembl", dataset = biomart_dataset)
entrezsymbol = biomaRt::getBM(attributes = c("entrezgene", "hgnc_symbol"), mart = mart)
entrezsymbol$entrezgene = as.character(entrezsymbol$entrezgene)
res = diff.exp.df
names(res) <- c("logFC", "AveExpr", "t", "P.Value", "adj.P.Val", "B", "hgnc_symbol")
enrich_rs = enrich_cp(res, "TNBC/LumA", type="all")
enrich_summary = enrich_rs$summary %>% arrange(p.adjust)
enrich_summary = convert_enriched_ids(enrich_summary,entrezsymbol = entrezsymbol) %>% arrange(p.adjust)
write.csv(enrich_summary, "datos/tnbc_luma_enrichment.csv")
Podemos ver que para los genes que se expresan diferencialmente, participan principalmente en el pathway de Neuroactive ligand-receptor interaction. Lo que es muy interesante, si buscamos en la literatura, podemos encontrar mucha información sobre este pathway.
El pathway neuroactive ligand-receptor interaction desempeña funciones en enfermedades neuronales, como por ejemplo en el caso del Parkinson. Aunque también se ha visto que desempeña funciones en el desarrollo de desórdenes adictivos, como es el caso de la dependencia al alcohol. Se han identificado cursos de tiempos en los que algunos genes pueden ser estimulados o inhibidos por estradiol en el cancer de mama humano, los resultados indican que estos cambios se enriquecen a las 24 horas de que un gen interacciones con este pathway. Un estudio de la subred en el cancer de mama reveló que los pathways oocyte meiosis eran subredes muy importantes de este pathway. Liu and Ye (n.d.)
Una vez hecho el enriquecimiento funcional de los genes, se realizó una selección de variables con el objetivo de quedarnos con los genes que resulten mas significativos. Para ello se ha dado la opción en el flujo de trabajo de:
# Sacamos mejores variables exhaustive
ranking.var<-function(df,indep.var){
indep.var<-"Surv(time,status)"
names.col<- colnames(df)[!colnames(df) %in% c("Node-Coded", "time", "status", "VitalStatus")]
pvalue<-c()
formula<-c()
for (i in names.col) {
ec.final<-as.formula(c("Surv(time,status) ~", i))
cox.data.norm<-coxph(ec.final, data = df)
cox.zph<-cox.zph(cox.data.norm)
#Comprobamos que la variable cumplen la proporcionalidad de riesgos
if(cox.zph$table[nrow(cox.zph$table),][3]>0.05){
formula<-c(formula, i)
pvalue<-c(pvalue,summary(cox.data.norm)$logtest["pvalue"])
}
}
variables<-ranking.var(clinicos.expresion, "Surv(time,status)")
formula<-as.simple.formula(paste(colnames(clinicos.expresion[!colnames(clinicos.expresion) %in% c("Node-Coded", "CompleteTCGAID","VitalStatus")]),collapse = "+"),"status")
weights <- random.forest.importance(formula, clinicos.expresion)
variables<-cutoff.k(weigths,50)
Me parecía muy interesante incluir aquí la selección de las variables de los modelos predictivos para poder hacer el análisis de enriquecimiento solo con los genes seleccionados. Si nos fijamos en los resultados del test de enriquecimiento para los datos reducidos.
Estos resultados podemos notar que no se parecen en nada a los anteriores, pero no dejan de ser interesantes. Si nos fijamos, podemos ver que los genes participan en la biosíntesis del ácido biliar. Esto resulta muy informativo, los pathways que asocian los ácidos biliares con el cancer suelen implicar un estrés oxidativo causando daños al DNA y generando inestabilidad genómica, apoptosis, …
Estos mecanismos también pueden ocurrir como efecto secundario de factores ambientales (estilo de vida, exposición a tóxicos,…) y su relación con cancer se ha identificado como crítica a diferentes niveles de tracto digestivo, así como a otros órganos como es el caso del cancer de mama. El efecto de cooridnado de diferentes agentes es posible (alcohol, fumar, …).Agostino Di Ciaula (n.d.)
Pero, cursiosamente también vemos que recientemente se ha descubierto que los ácidos biliares también pueden tener efectos contra el cancer, ya que modulan los mismos pathways que la toxicidad, la apoptosis, …
Antes de replicar los datos, vamos a montar el data frame con los genes que hemos seleccionado en el análsis de enriquecimiento. Para replicar los datos he creado una función muy sencilla en la que podemos establecer el step. Y nos devuelve los datos replicados.
replicate.cum <- function(step, data, min.time, max.time){
data.cum <- data.frame()
pb<-txtProgressBar(min=0, max =dim(clinicos.expresion)[1], style = 3 )
for( i in 1:dim(clinicos.expresion)[1]){
linea <- data[i,]
for( j in seq(min.time, max.time, by=step)){
lin <- linea
if(as.numeric(linea$`OS Time`) >j){
lin$`OS Time` <- j
lin$`OS event` <- 0
}else{
if(as.numeric(linea$`OS event`)==2){
lin$`OS Time` <- j
lin$`OS event` <- 1
}
}
data.cum <- rbind(data.cum, lin)
}
setTxtProgressBar(pb,i)
}
return(data.cum)
}
rep.data <- replicate.cum(step = step, clinicos.expresion,
min.time = as.numeric(min(clinicos.expresion$`OS Time`, na.rm = TRUE)),
max.time = as.numeric(max(clinicos.expresion$`OS Time`, na.rm = TRUE)))
Esta forma de replicar los datos pese a que puede parecer normal, provoca ciertos fallos. Concretamente, me genera un dataset con los datos desbalanceados. Es decir, que para status tengo muchos mas 0 que 1.
Al tratar de ejecutar modelos de machine learning con datos desbalanceados, como me pasó a mi, lo que ocurre es que el modelo a la hora de predecir me devuelve siempre la misma clase, y como en el test también tenemos muchos mas 0 que 1 puede parecer que estamos teniendo buenos resultados. Haciendo pruebas para ver que estaba pasando, modificando datos y demás obtenía resultados como este:
kable(as.data.frame(cm$byClass)) %>%
kable_styling()
| cm$byClass | |
|---|---|
| Sensitivity | 1.0000000 |
| Specificity | 0.2000000 |
| Pos Pred Value | 0.9875260 |
| Neg Pred Value | 1.0000000 |
| Precision | 0.9875260 |
| Recall | 1.0000000 |
| F1 | 0.9937238 |
| Prevalence | 0.9844560 |
| Detection Rate | 0.9844560 |
| Detection Prevalence | 0.9968912 |
| Balanced Accuracy | 0.6000000 |
Pero si nos fijamos detenidamente obtenemos valores para el Balanced Accuracy tenemos un 0.6.
Hay varios métodos de muestreo que han sido diseñados para tratar datos desbalanceados, los cuales pueden ser agrupados en cuatro categorías:
Todos estos métodos modifican la proporción de las clases y el tamaño del dataset original. Yo he decidido utilizar un método de submuestreo. Para ello, eliminaremos observaciones de la clase mayoritaria con el fin de igualar los tamaños de las clases.
Apliquemos la función undersample_ds al dataset training requiriendo un número de 400 observaciones para cada clase.
undersample_ds <- function(x, classCol, nsamples_class){
for (i in 1:length(unique(x[, classCol]))){
class.i <- unique(x[, classCol])[i]
if((sum(x[, classCol] == class.i) - nsamples_class) != 0){
x <- x[-sample(which(x[, classCol] == class.i),
sum(x[, classCol] == class.i) - nsamples_class), ]
}
}
return(x)
}
training_bc <- undersample_ds(rep.data, "status", nsamples_class)
table(rep.data.balanced$status)
##
## 0 1
## 623 623
Teniendo los datos balanceados obtenemos los siguientes resultados:
kable(as.data.frame(cm1$byClass)) %>%
kable_styling()
| cm1$byClass | |
|---|---|
| Sensitivity | 0.9850746 |
| Specificity | 1.0000000 |
| Pos Pred Value | 1.0000000 |
| Neg Pred Value | 0.9830508 |
| Precision | 1.0000000 |
| Recall | 0.9850746 |
| F1 | 0.9924812 |
| Prevalence | 0.5360000 |
| Detection Rate | 0.5280000 |
| Detection Prevalence | 0.5280000 |
| Balanced Accuracy | 0.9925373 |
Podemos ver que el accuracy balanceado ahora es muy similar al accuracy normal. Estos resultados que se han mostrado son para el mejor alpha obtenido entrenando una elasticnet con todas las variables.
Para esta parte se ha empezado con una doble validación cruzada que hemos implementado a mano basándonos en la función que ofrece el paquete glmnet cv.glmnet.
foldDoubleKVal<-function(df){
all.train.test <- sparse.model.matrix(as.formula(paste("status~", paste(colnames(df)[colnames(df)!="status"], sep="", collapse = "+"))),df)
pb<-txtProgressBar(min=1, max =10, style = 3 )
# Input 2. Set the fractions of the dataframe you want to split into training,
# validation, and test.
fractionTraining <- 0.60
fractionValidation <- 0.20
fractionTest <- 0.20
# Compute sample sizes.
sampleSizeTraining <- floor(fractionTraining * nrow(all.train.test))
sampleSizeValidation <- floor(fractionValidation * nrow(all.train.test))
sampleSizeTest <- floor(fractionTest * nrow(all.train.test))
auc.list<-c()
accuracy<-c()
accuracy.balanced<-c()
best.alpha<-c()
best.lambda<-c()
for(j in 1:10) {
# Create the randomly-sampled indices for the dataframe. Use setdiff() to
# avoid overlapping subsets of indices.
indicesTraining <- sort(sample(seq_len(nrow(all.train.test)), size=sampleSizeTraining))
indicesNotTraining <- setdiff(seq_len(nrow(all.train.test)), indicesTraining)
indicesValidation <- sort(sample(indicesNotTraining, size=sampleSizeValidation))
indicesTest <- setdiff(indicesNotTraining, indicesValidation)
# Finally, output the three dataframes for training, validation and test.
x.train <- all.train.test[indicesTraining, ]
x.val <- all.train.test[indicesValidation, ]
x.test <- all.train.test[indicesTest, ]
y.train <- df[indicesTraining, ]$status
y.val <- df[indicesValidation, ]$status
y.test <- df[indicesTest, ]$status
# Vemos el alpha para el que obtenemos el mejor resultado
lista<-list()
aucs<-c()
lambdas<-c()
for (i in 0:10) {
lambdas<-c(lambdas,cv.glmnet(x.train, y.train, family = "binomial",
nfold = 10, type.measure = "auc", paralle = TRUE, alpha = 1)$lambda.1se)
lista[[(i+1)]]<-glmnet(x.train, y.train, family = "binomial",lambda = lambdas[i+1],alpha = i/10)
aucs<-c(aucs,roc(y.val, as.numeric(predict(lista[[(i+1)]], x.val, type = "response"))))
}
best.alpha<-c(best.alpha, (which.min(lambdas)-1)/10)
best.lambda<-c(best.lambda, min(lambdas))
lasso.model<-glmnet(x.train, y.train, family = "binomial",lambda = min(lambdas),alpha = (which.min(lambdas)-1)/10)
lasso.prob <- predict(lasso.model,type="class", newx = x.test)
#Sacamos la matriz de confusion
cm = confusionMatrix(as.factor(lasso.prob),as.factor(y.test))
accuracy.balanced<-c(accuracy.balanced,cm$byClass["Balanced Accuracy"])
accuracy<-c(accuracy, cm$byClass["Balanced Accuracy"])
auc.list<-c(auc.list, roc(y.test, as.numeric(lasso.prob))$auc[1])
setTxtProgressBar(pb,j)
}
return(data.frame(auc=auc.list, accuracy= accuracy, balanced.accuracy = accuracy.balanced, alpha = best.alpha, lambda = best.lambda))
}
resDoubleValidation<-foldDoubleKVal(df)
Los resultados los hemos guardado usando la función kable.
loadpkg("kableExtra")
loadpkg("magick")
kable(resDoubleValidation, format = "latex")%>%
kable_styling() %>%
save_kable(file="results/Lasso-ridge/10DoubleValidation.jpg")
cat("\n The mean AUC for lasso/ridge is ", mean(resDoubleValidation$accuracy))
##
## The mean AUC for lasso/ridge is 0.9584192
En esta parte hemos guardado los resultados para después compararlos con los obtenidos con el resto de modelos de machine learning.
Además he añadido unas gráficas sobre el entrenamiento de la elasticnet para 1=LASSO, 0=RIDGE y 0.5, que están guardados como todos los resultados en la carpeta results/Lasso-ridge.
Pero nos llama la atención los resultados tan sorprendentemente buenos que arroja el modelo. De no ser por un error, estos resultados solo se podrían explicar por que tenemos una fuerte dependencia lineal entre los datos de entrenamiento y los resultados. Me imaginaba que al balancear los datos obtendría los valores de accuracy que estaba obteniendo antes para el accuracy balanceado.
Por lo que voy a utilizar 3 modelos de machine learning para discutir los resultados:
Se ha implementado una función para calcular el 10 Cross Validation AUC.
crossValidation.models <-function(data,k) {
auc.glm<-list()
auc.dt<-list()
auc.nnet<-list()
cf<-createFolds(data$status, k=k)
formula=as.formula(paste("status", ".", sep="~"))
for(i in 1:k) {
test <-data[unlist(cf[i]),]
train <-data[-unlist(cf[i]),]
train$status <-as.factor(train$status)
test$status <-as.factor(test$status)
(l <- sapply(train, function(x) is.factor(x)))
m <- train[, l]
ifelse(n <- sapply(m, function(x) length(levels(x))) == 1, "DROP", "NODROP")
#dt
modelo.dt<-rpart(formula,data=train,control=rpart.control(minsplit=1, cp=0.007));
mod.estimado.dt<-predict(modelo.dt, newdata=test, type="matrix");
mod.estimado.dt<-as.matrix(mod.estimado.dt)
if(is.factor(mod.estimado.dt)) mod.estimado.dt<-ifelse(mod.estimado.dt=="X1", 1, 0)
new.auc.dt<-roc(as.numeric(test$status), as.numeric(mod.estimado.dt[,dim(mod.estimado.dt)[2]]), smooth=F, auc=T)$auc
auc.dt<-c(auc.dt,new.auc.dt)
#Regresión logística
modelo.glm<-glm(formula, train, family=binomial("logit"));
modelo.glm$xlevels$PAM50mRNA <- union(modelo.glm$xlevels$PAM50mRNA, levels(test$PAM50mRNA))
mod.estimado.glm<-predict(modelo.glm, newdata=test, type="response");
if(is.factor(mod.estimado.glm)) mod.estimado.dt<-ifelse(mod.estimado.glm=="X1", 1, 0)
new.auc.glm<-roc(as.numeric(test$status), mod.estimado.glm, smooth=F, auc=T)$auc
auc.glm<-c(auc.glm,new.auc.glm)
#Neural network
modelo.nnet<-nnet(formula,data=train, size=5, decay=0, maxit=7, trace=F);
mod.estimado.nnet<-predict(modelo.nnet, newdata=test, type="raw");
if(is.factor(mod.estimado.nnet)) mod.estimado.dt<-ifelse(mod.estimado.nnet=="X1", 1, 0)
new.auc.nnet<-roc(as.numeric(test$status), mod.estimado.nnet, smooth=F, auc=T)$auc
auc.nnet<-c(auc.nnet,new.auc.nnet)
}
all.means<-data.frame(method=c("glm","dt","nnet"), metrics=c(mean(unlist(auc.glm)), mean(unlist(auc.dt)), mean(unlist(auc.nnet))))
return(all.means)
}
debug(crossValidation.models)
res<-crossValidation.models(rep.data.balanced,8)
Se ha decidido no añadir modelos de Deep Learning porque los resultados para los modelos mas sencillos están dando mejores resultados.
Para esta parte vamos a hacer 2 veces survival, una con todos los pacientes para ver la supervivencia global y otra con los pacientes clusterizados. Y compararemos los resultados. Para el flujo de trabajo, simplemente se van a guardar las diferentes curvas de supervivencia.
#Estratificamos las variables necesarias
df$age<-cut(as.numeric(df$`AgeatInitialPathologicDiagnosis`), breaks = c(0,35,60,100))
df$time<-as.numeric(df$time)
df$status<-as.numeric(df$status)
df$AJCC<-as.factor(df$AJCCStage)
# Vemos las estratificaciones de las variables de interes
fit.age<-survfit(Surv(time, status) ~ age, data=df)
fit.gen<-survfit(Surv(time, status) ~ Gender, data=df)
fit.node<-survfit(Surv(time, status) ~ Node, data=df)
fit.meta<-survfit(Surv(time, status) ~ Metastasis, data=df)
fit.ajcc<-survfit(Surv(time, status) ~ AJCC, data=df)
#Ploteamos los resultados
splots<-list()
splots[[1]] <- ggsurvplot(fit.age, censor= TRUE, main="Estratificado: edad", pval = TRUE)
splots[[2]] <- ggsurvplot(fit.gen, censor= TRUE, main="Estratificado: genero", pval = TRUE)
splots[[3]] <- ggsurvplot(fit.node, censor= TRUE, main="Estratificado: nodos", pval = TRUE)
splots[[4]] <- ggsurvplot(fit.meta, conf.int = TRUE, censor= TRUE, main="Estratificado: Metastasis", pval = TRUE)
splots[[5]] <- ggsurvplot(fit.ajcc, censor= TRUE, main="Estratificado: Metastasis", pval = TRUE)
En el codigo superior, lo que hacemos es estratificar una serie de variables que pueden resultar de interes, y mostramos la curva de supervivencia estratificando para cada uno de los niveles seleccionados en cada variable.
En las estratificaciones globales podemos observar como el mejor resultado se obtiene con diferencia para la variable de nodos. Vemos que en función del AJCC se generan curvas muy diferentes entre si, pero que no aparecen diferenciadas significativamente. Y no encontramos tampoco diferencias significativa entre las edades, solo se aprecia que desde 35 hacia adelante tenemos un grupo de pacientes pero para pacientes con edades inferiores a 35 no tenemos suficientes datos.
Hierarchical clustering es una alternativa a los métodos de partitioning clustering que no requiere que se pre-especifique el número de clusters. Los métodos que engloba el hierarchical clustering se subdividen en dos tipos dependiendo de la estrategia seguida para crear los grupos:
Agglomerative clustering (bottom-up): el agrupamiento se inicia en la base del árbol, donde cada observación forma un cluster individual. Los clusters se van combinado a medida que la estructura crece hasta converger en una única “rama” central.
Divisive clustering (top-down): es la estrategia opuesta al agglomerative clustering, se inicia con todas las observaciones contenidas en un mismo cluster y se suceden divisiones hasta que cada observación forma un cluster individual.
En ambos casos, los resultados pueden representarse de forma muy intuitiva en una estructura de árbol. Los métodos de clustering son unsupervised, lo que significa que al aplicarlos no se hace uso de la variable respuesta. Una vez obtenidos los resultados se añade esta información para determinar si es posible agrupar a las líneas celulares empleando su perfil de expresión genico.
expresion <- scale(df, center = TRUE, scale = TRUE)
matriz_distancias <- dist(x = expresion, method = "euclidean")
hc_completo <- hclust(d = matriz_distancias, method = "complete")
hc_average <- hclust(d = matriz_distancias, method = "average")
hc_single <- hclust(d = matriz_distancias, method = "single")
La elección del tipo de linkage influye de forma notable en los dendrogramas obtenidos. Por lo general, single linkage tiende a formar clusters muy grandes en los que las observaciones individuales se unen una a una. El completo y average suele generar dendrogramas más compensados con clusters más definidos, tal como ocurre en este ejemplo.
En la imagen superior tenemos los silhouettes para diferentes K’s. Tras ver varios resultados he decidido coger un numero de clusters igual a 2.
clusters <- cutree(tree = hc_completo, k = 2)
La variable clusters ahora la vamos a utilizar para dividir el set de datos y volver a revisar las curvas de supervivencia.
group.1<-df[clusters,]
group.2<-df[-clusters,]
fit.age.1<-survfit(Surv(time, status) ~ age, data=group.1)
fit.node.1<-survfit(Surv(time, status) ~ Node, data=group.1)
fit.ajcc.1<-survfit(Surv(time, status) ~ AJCC, data=group.1)
fit.age.2<-survfit(Surv(time, status) ~ age, data=group.2)
fit.node.2<-survfit(Surv(time, status) ~ Node, data=group.2)
fit.ajcc.2<-survfit(Surv(time, status) ~ AJCC, data=group.2)
Para las variables que hemos visto de interes anteriormente, vamos a volver a hacer el análisis de supervivencia.
# Vemos las estratificaciones de las variables de interes
fit.age.1<-survfit(Surv(time, status) ~ age, data=group.1)
fit.node.1<-survfit(Surv(time, status) ~ Node, data=group.1)
fit.ajcc.1<-survfit(Surv(time, status) ~ AJCC, data=group.1)
fit.age.2<-survfit(Surv(time, status) ~ age, data=group.2)
fit.node.2<-survfit(Surv(time, status) ~ Node, data=group.2)
fit.ajcc.2<-survfit(Surv(time, status) ~ AJCC, data=group.2)
splots1<-list()
splots1[[1]] <- ggsurvplot(fit.age.1, censor= TRUE, main="Estratificado: edad", pval = TRUE)
splots1[[2]] <- ggsurvplot(fit.node.1, censor= TRUE, main="Estratificado: nodos", pval = TRUE)
splots1[[3]] <- ggsurvplot(fit.ajcc.1, censor= TRUE, main="Estratificado: AJCC", pval = TRUE)
splots1[[4]] <- ggsurvplot(fit.age.2, censor= TRUE, main="Estratificado: edad", pval = TRUE)
splots1[[5]] <- ggsurvplot(fit.node.2, censor= TRUE, main="Estratificado: nodos", pval = TRUE)
splots1[[6]] <- ggsurvplot(fit.ajcc.2, censor= TRUE, main="Estratificado: AJCC", pval = TRUE)
Podemos notar claramente que el cluster se encuentra desbalanceado, tenemos muchos mas pacientes a un lado que en otro y esto hace que un cluster se parezca mucho a la original, mientras que el otro no tenga datos suficientes como para poder compararlo.
Me parece importante comentar los resultados que he ido teniendo a lo largo de la experimentación. Dejando de lado los test de enriquecimiento, que creo que están bien documentados, la verdad es que los resultados que se han ido obteniendo parecen cuanto menos extraños.
Agostino Di Ciaula, † Emilio Molina-Molina, * David Q.-H. Wang. n.d. Bile Acids and Cancer:Direct and Environmental-Dependent Effects. "http://www.medigraphic.com/pdfs/hepato/ah-2017/ahs171k.pdf".
Liu, Huijuan, and Hui Ye. n.d. Screening of the Prognostic Targets for Breast Cancer Based Co-Expression Modules Analysis. "https://www.ncbi.nlm.nih.gov/pmc/articles/PMC5646985/".